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Abstract 

The  focus  of  this  research  was  to  derive  a  new  algorithm  for  correction  of  gain 
nonuniformities  in  EIDAR  focal  plane  arrays  using  as  few  frames  as  possible.  Because 
of  the  current  low  production  rate  of  EIDAR  focal  plane  arrays  there  is  a  natural  tendency 
for  extreme  nonuniformities  to  exist  on  a  pixel  by  pixel  basis  as  the  manufacturing 
technique  has  not  yet  been  perfected.  Generally,  nonuniformity  correction  techniques 
require  a  large  number  of  frames  and/or  have  obscure  requirements  on  the  translational 
shifts  in  the  input  image  frames.  This  thesis  presents  a  solution  for  finding  multiplicative 
nonuniformities  that  exist  in  a  focal  plane  array  and  mitigating  the  effect  those 
nonuniformities . 
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FAST  SCENE  BASED  NONUNIFORMITY  CORRECTION  WITH  MINIMAL 


TEMPORAL  LATENCY 


I.  Introduction 


LIDAR  (Light  Detection  and  Ranging),  which  uses  the  same  basic  principles  of 
RADAR  (Radio  Detection  and  Ranging),  relies  on  a  pulse  of  light  and  a  sensor  usually 
consisting  of  a  EPA  (Local  Plane  Array)  consisting  of  a  set  of  detectors.  This  EPA  views 
a  scene  of  reflected  light  from  objects  in  the  scene.  These  light  returns  from  objects  are 
then  quickly  processed  to  determine  the  range  these  pulses  were  returned  from  based  on 
the  LIDAR  system  position.  This  thesis  specifically  focuses  on  LIDAR  EPA  sensors  as 
the  algorithm  requires  a  scene  with  at  least  a  full  pixel  of  motion  between  frames  to 
mitigating  the  effect  of  the  spatial  nonuniformities  in  the  sensor  EPA.  A  LIDAR  example 
is  shown  in  figure  1  below: 


Eigure  1  -  An  eagle  eye  view  of  a  LIDAR  system  (I)  imaging  three  objects  (A,  B,  and  C) 
and  a  coherent  pulse  for  three  different  times. 


As  mentioned  before,  the  FPAs  for  any  type  of  imager  eonsist  of  an  array  of 
deteetors,  eaeh  deteetor  of  whieh  ean  also  be  referred  to  as  a  pixel.  Eaeh  deteetor 
responds  to  photons  that  eollide  with  it  by  generating  a  eharge.  This  eharge  is  then 
amplified  and  eonverted  from  a  voltage  to  a  digital  signal  using  a  readout  integrated 
eireuit  (ROIC)  resulting  in  numbers  of  digital  units  (DU)  proportional  to  the  intensity  of 
the  light  ineident  on  the  deteetor.  This  is  done  for  eaeh  deteetor  in  the  FPA. 

In  new  and  old  designs  of  FPAs  alike,  when  an  FPA  is  manufaetured  there  are 
ineonsisteneies  between  the  gain  and  bias  of  eaeh  pixel  in  the  eonversion  from  photons  to 
a  DU.  Note  that  these  ineonsisteneies  also  ehange  as  time  and/or  sensor  temperature 
ehange,  giving  rise  to  the  requirement  of  frequent  ealibration  in  FPAs.  This  is  espeeially 
true  for  FPAs  of  newer  design,  as  the  manufaeturing  proeess  for  newer  types  of  FPAs  has 
yet  to  be  perfeeted.  This  is  eurrently  true  for  FIDAR  FPAs  as  their  design  and 
manufaeturing  proeess  is  still  relatively  new  at  the  time  this  was  presented. 

For  illustration  purposes,  imagine  the  unlikely  event  that  the  same  number  of 
photons  are  ineident  on  eaeh  deteetor.  Beeause  of  the  nature  of  these  deteetors,  eaeh  one 
will  generate  a  different  level  of  eharge.  If  these  nonuniformities  in  pixel  response  do  not 
ehange  over  time  eaeh  deteetor  eould  be  ealibrated  onee  after  the  manufaeturing  proeess 
and  used  indefinitely  without  error.  Unfortunately  this  is  not  the  ease.  Over  time  these 
nonuniformities  in  the  FPA  ehange  and  the  images  returned  from  the  FPA  need  to  be 
ealibrated  frequently. 

Calibration  to  remove  nonuniformities,  regardless  of  sensor  type,  ean  be  a  very 
expensive  proeess  and  may  also  result  in  data  being  missed  by  the  sensor  beeause  a 
ealibration  eyele  needed  to  take  plaee.  An  ideal  ealibration  would  be  one  that 
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continuously  calibrates  the  sensor  array  based  on  typieal  imagery  input  so  that  no  desired 
imagery  would  be  missed  and  no  down  time  for  the  sensor  would  be  required. 

As  a  general  rule,  visible,  infrared,  and  LIDAR  sensor  arrays  alike  are  all 
suseeptible  to  nonuniformities  in  varying  degrees  that  ehange  over  time  and  all  generally 
need  some  type  of  eorreetion. 


The  Need  for  LIDAR  Nonuniformity  Correction 

As  mentioned  previously,  3-D  LIDAR  FPAs  are  a  relatively  new  sensor  paekage 
and  are  eurrently  in  very  low  produetion  eompared  to  other  imager  types  available. 
Beeause  of  this,  the  manufaeturing  proeess  for  these  sensor  arrays  is  still  being  perfeeted, 
yielding  arrays  with  drastic  nonuniformities  notieeable  to  even  the  naked  eye.  Beeause  of 
the  degree  of  these  nonuniformities,  there  is  an  immediate  need  for  an  algorithm  that  will 
reduee  nonuniformities  in  eurrent  LIDAR  FPAs. 

3-D  LIDAR  sensors,  sometimes  referred  to  as  3D  eameras  or  flash  LIDAR,  are 
generally  cameras  with  an  FPA  that  ean  proeess  an  unusually  high  number  of  frames  over 
a  very  short  period  of  time.  For  example,  the  “Zenith  eamera”  refereneed  in  Elkins  et  al. 
ean  proeess  16  frames  at  a  speed  of  100  million  frames  per  second  (Mfps)  [1].  This 
sensor  has  a  very  quiek  aequisition  cyele,  but  a  very  slow  data  readout  eyele,  and  is 
referred  to  as  a  fast-in-slow-out  (FISO)  proeess.  This  is  currently  true  for  nearly  all 
eurrent  LIDAR  arrays  due  to  the  vast  amounts  of  data  yielded  from  eolleeting  range  data 


for  eaeh  deteetor  in  the  FPA. 


The  potential  future  use  for  LIDAR  sensors  ean  be  been  seen  in  reeent  finaneial 
investments  and  published  artieles  [6]  [11].  Flash  LIDAR  sensors  have  many  possible 
future  applieations,  and  all  would  benefit  from  the  removal  of  nonuniformities  from  the 
resulting  data.  Additionally,  beeause  of  the  nature  of  LIDAR  FPAs  it  is  eurrently 
difficult  and  expensive  to  retrieve  flat  field  data  in  order  to  perform  a  complete 
calibration.  Because  of  this,  an  ideal  algorithm  could  calibrate  the  LIDAR  FPA  with  data 
that  needs  to  be  collected  already,  or  is  easily  added  to  a  collection  schedule,  reducing 
overall  costs. 

The  algorithm  described  in  this  document  is  designed  to  reduce  gain 
nonuniformities  while  only  requiring  a  small  number  of  frames,  thereby  giving  the  ability 
to  use  data  that  may  have  already  been  collected.  Because  the  scene  can  be  captured 
without  returns  from  the  LIDAR  pulse,  the  bias  can  be  found  for  each  scene  observed. 
With  the  bias  taken  care  of,  the  algorithm  presented  here  can  focus  on  the  sensor  gain 
values  for  each  detector.  This  makes  the  algorithm  especially  useful  and  applicable  to 
LIDAR  data. 


II,  Uniqueness  of  LIDAR  Nonuniformity  Correction  Algorithm 

There  are  many  technical  articles,  old  and  new,  discussing  the  methods  for 
removing  nonuniformities  from  the  resulting  image  data  collected  by  FPAs  of  all  types, 
including  visible  and  infrared.  However,  a  comprehensive  literature  search  suggests 
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there  are  very  few,  if  any,  articles  that  discuss  the  removal  of  nonuniformities  from 


LIDAR  FPAs. 


Literature  Search  for  LIDAR  Nonuniformity  Correction  Techniques 

A  search  through  journals  and  conference  proceedings  show  that  there  are 
presently  no  articles  that  discuss  methods  for  correction  of  nonuniformities  in  LIDAR 
FPAs.  There  are  many  articles  that  discuss  methods  for  the  removal  of  nonuniformities 
in  other  types  of  FPAs  as  referenced  ([4],  [9],  and  [3]),  many  specifically  concerned  with 
infrared  FPAs,  however  there  is  some  discontinuity  when  trying  to  apply  these  types  of 
algorithms  to  image  data  from  a  LIDAR  collection. 

First,  LIDAR  image  data  comes  in  cubes,  which  can  be  compared  to  video  data 
from  an  infrared  or  visible  imager.  However,  the  primary  difference  is  that  the  data  in  a 
LIDAR  cube  has  a  piece  of  a  scene  in  each  frame  corresponding  to  a  range  bin  for  an 
equal  distance  in  the  scene,  whereas  each  frame  in  a  data  cube  from  a  standard  imager 
contains  intensities  of  the  entire  scene  integrated  over  a  longer  time,  or  over  many  range 
bins,  equivalent  to  the  depth  of  the  scene,  removing  distance  information  from  data 
collection.  This  fact  must  be  kept  in  mind  when  constructing  an  algorithm  or  evaluating 
algorithms  that  are  used  to  remove  nonuniformities  from  typical  FPAs. 

Second,  most  algorithms  in  existence  remove  nonuniformities  from  FPAs  using  a 
large  number  of  frames  to  iteratively  solve  for  gain  and  bias.  In  a  LIDAR  sensor  the  bias 
is  already  recorded  in  the  frames  before  and  after  the  return  of  the  illuminating  pulse  for 
the  entire  scene,  with  some  additional  background.  This  provides  an  opportunity  that 
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would  not  normally  be  available  for  a  standard  visible  or  infrared  imager  beeause  there  is 
no  absenee  of  the  illuminated  light  souree  as  there  is  in  a  LIDAR  sensor  system. 


Comparison  to  Current  Methods 

There  are  primarily  two  methods  for  nonuniformity  eorreetion,  seene-based 
eorreetion  and  referenee-based  eorreetion  [4], 

Referenee-based  algorithms  use,  as  implied,  referenees  to  calibrate  the  sensors  in 
the  FPA  [9].  This  method  of  calibration  uses  multiple  collections  consisting  of  differing 
calibration  targets  to  develop  a  model  for  sensor  response,  usually  linear.  While  the 
algorithm  discussed  in  this  thesis  is  somewhat  similar  to  the  algorithm  discussed  in  [9], 
reference-based  algorithms  are  too  dissimilar  to  allow  a  meaningful  comparison  because 
they  only  use  multiple  calibrated  scenes  in  order  to  calibrate  the  FPA  and  this  LIDAR 
based  algorithm  uses  the  bias  reference  at  the  beginning  or  end  of  a  LIDAR  image 
sequence  to  aid  in  calibration.  The  algorithm  described  here  does  not  require  calibrated 
targets  but  does  rely  on  scene  shifts  (Scene  shifts  being  shifts  of  the  viewed  scene  in  the 
vertical  and/or  horizontal  direction.)  during  collection  and  the  resulting  image  data 
including  those  shifts. 

Scene-based  algorithms,  or  algorithms  that  rely  solely  on  scene  data  collected  by 
the  imager,  are  only  slightly  more  comparable  to  the  algorithm  described  in  this  paper. 
Specifically,  only  scene-based  algorithms  that  require  a  low  number  of  frames  could  be 
compared  to  this  algorithm  because  algorithms  that  require  high  numbers  of  frames  are 
not  applicable  to  LIDAR  data  containing  a  relatively  small  number  of  observations. 
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Additionally,  some  current  scene-based  algorithms  tend  to  require  an  image  sequence  that 
exhibits  a  restrained  set  of  shifts  required  by  the  algorithm  and  a  satisfactory  result  can 
not  be  achieved  unless  the  input  frames  satisfy  these  restraints.  For  example,  one 
algorithm  may  require  image  pairs  that  exhibit  a  one  dimensional  shift  in  both  the 
horizontal  and  the  vertical  direction  individually  and  image  pairs  that  exhibit  shifts  in 
both  the  horizontal  and  the  vertical  direction  simultaneously  for  the  algorithm  to  run  with 
satisfactory  results,  as  referenced  by  Hayat  et  al  [3].  The  algorithm  described  by  Hayat  et 
al.  only  corrects  bias  and  not  gain  whereas  the  algorithm  presented  here  addresses  the 
gain  and  not  the  bias. 

These  differences,  and  the  fact  that  no  current  nonuniformity  correction 
algorithms  have  been  demonstrated  with  LIDAR  FPAs,  make  this  algorithm  new  and 
difficult  to  compare  to  existing  algorithms. 
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Ill,  Derivation  of  Nonuniformity  Correction  Algorithm 


In  order  to  derive  a  nonuniformity  eorreetion  algorithm  that  is  not  dependent  on  a 
high  number  of  frames,  eaeh  pixel  is  modeled  with  an  estimated  linear  response.  This  is 
generally  a  very  reasonable  approximation  and  eommon  when  estimating  a  response  from 
a  deteetor. 

Also,  this  algorithm  solves  for  the  gain,  not  the  bias,  as  all  seene  based 
nonuniformity  eorreetion  algorithms  that  use  a  small  number  of  frames  do,  beeause  the 
bias  in  a  LIDAR  sensor  ean  be  seen  in  the  eolleetion  of  data  just  before  responses  eome 
baek  from  the  emitted  pulse  in  a  eolleetion  sequenee,  or  just  thereafter.  Using  this 
information  that  would  not  normally  be  available  to  a  standard  FPA,  a  formulation  of  the 
gain  is  iterated  to  a  solution  in  order  to  remove  the  nonuniformities  in  the  FPA. 


Derivation  of  Algorithm 

Assumptions  must  be  made  in  order  to  derive  an  algorithm  to  eorreet  the 
nonuniformities  in  a  LIDAR  array.  It  is  assumed  that  an  image  is  gathered  at  a  time 
when  no  pulse  is  present.  It  is  also  assumed  that  the  response  for  eaeh  deteetor  is  linear. 

Additionally,  beeause  of  the  natural  Poisson  nature  of  photon  arrivals  [2],  the 
PMF  (Probability  Mass  Funetion)  of  the  natural  proeess  of  eounting  photons  can  be 
described  generally  with  a  Poisson  random  variable,  n,  as  seen  in  Leon-Garcia  [5]; 
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Where  in  equation  (3.1)  n  is  the  mean  number  of  photons  per  measurement  interval  and 
is  the  PMF.  Note  that  in  general  the  log-likelihood,  or  natural  log  of  the  PMF  for  a 
single  measurement,  can  be  described  with  the  following  equation: 

L{n,n)  =  n-\n{n)-n-\n{n)  (3.2) 

In  order  to  maximize  the  log-likelihood  of  a  function,  the  derivative  of  the  log- 
likelihood  is  taken  with  respect  to  n  and  set  equal  to  zero,  knowing  that  a  maximization 
of  the  log  of  a  function  will  maximize  the  function  as  well: 

=  (3.3) 

d{n)  n 

For  further  information,  these  estimation  techniques  are  described  in  Van  Trees 
starting  on  page  63  [10],  including  a  brief  description  of  the  log-likelihood. 

With  these  ideas  in  mind,  the  algorithm  is  derived  by  first  formulating  the  model 
of  the  data: 

d (x,  y)  =  a{x,  y)  ■  z(x,  y)  +  n(x,  y)  (3 .4) 

Where  in  equation  (3.4)  d{x,y)  is  the  resulting  image  with  nonuniformities  along 
the  two  dimensions  X,  y  and  are  both  integers,  a(x,y)  are  the  nonuniformities  along  the 
array,  i{x,y)  is  the  original  scene  the  algorithm  is  attempting  to  recreate,  and  n{x,y)  is 
some  noise  associated  with  this  image  frame  that  makes  d(x,y)  Poisson  in  nature.  Next, 
in  order  to  account  for  shifted  data  the  shifts  must  be  included  in  the  image  model: 

dk  {x,  y)  =  a{x,  y)-i{x- j3^,y-s^)  +  n^  (x,  y)  (3.5) 

Where  in  equation  (3.5)  dj^(x,y)  is  the  resulting  image  for  the  k*  frame, 
i(x-  J3^,y-s^)  is  the  original  scene  shifted  inx  by  /d^.  andy  by  s,^  where  all  shifts  are 
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integers,  and  nj^{x,y)  is  the  noise  associated  with  the  A:*  frame.  Using  this  notation,  two 
different  frames  could  be  described  by  the  following  set  of  equations: 

dy{x,y)  =  a{x,y)-i{x- py,y  -  s^)  +  ny{x,y) 

(3.6) 

d2{x,y)  =  a{x,y)-i{x-fi2,y-^2)  +  n^{x,y) 

Also,  a  shifted  image  can  be  described  as  the  convolution  sum  over  dummy 
variables  (y,//)  between  the  image  and  a  shifted  Dirac  delta  function  as  shown  in 
equation  (3.7),  this  will  be  used  later: 

Kx- P„y-£,)  =  Yj'Yj  Kr,'n)-d{x-r- p„y-ri-s,)  (3.7) 

7  1 

Given  this  notation,  the  iterative  result  is  found  by  restraining  the  sum  of  i  to  be 
one.  This  is  accomplished  via  a  Lagrange  multiplier.  A, .  This  A,  multiplies  the  sum  of  i 
and  is  subtracted  from  the  log-likelihood  function.  An  example  of  the  Lagrange  being 
used  is  given  in  Van  Trees  starting  on  page  33  [10].  The  constraint  is  then  that  the 
average  of  the  image  must  be  equal  to  the  average  of  the  data: 

'  - =  (3-8) 

X  y 

Where  K  is  the  total  number  of  frames  used.  Given  that  the  scene  itself  can  be  described 
as  a  random  variable  with  a  joint  distribution  of  each  of  the  Poisson  pixel  PMFs,  the  PMF 
of  such  a  random  variable  can  be  described  as  follows  in  terms  of  a  pixel  in  x  and  y: 

P(D(x,y)  =  d(x,y))  =  e-''"  (3.9) 

y)  • 
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Given  this  PMF  for  a  pixel,  the  joint  PMF  for  the  image  ean  be  deseribed  by  the 
produet  of  the  individual  pixel  PMFs  if  all  the  measurements  at  eaeh  pixel  are  statistieally 
independent; 


(3.10) 


In  order  to  remove  the  nonuniformities,  the  log-likelihood  funetion  of  a  Poisson 
PMF  must  take  the  form  described  above  in  equation  (3.10).  Because  the  natural 
logarithm  is  a  monatomic  function,  maximizing  the  log-likelihood  will  also  maximize  the 
joint  PMF  of  the  data.  A  maximum  likelihood  (ML)  estimate  can  be  obtained  in  order  to 
remove  the  nonuniformities  from  the  original  data.  Note  that  ML  is  an  acceptable 
technique  in  this  scenario  because  there  is  no  knowledge  of  the  statistics  of  the 
nonuniformities  and  therefore  a  maximum  a  posteriori  (MAP)  estimate  cannot  be  used. 
Additionally,  minimum-mean  square  error  (MMSE)  estimation  cannot  be  used  because 
the  priors  are  not  known,  making  ML  a  good  overall  choice  as  an  estimation  scheme. 
Using  the  form  for  the  log-likelihood  above,  the  log-likelihood  for  the  image  model  is 
formed  for  all  the  pixels  in  the  image; 

L(a,i)  =  III  (x,  y)  ■  ln[«(x,  y)  •  i(x  -  ,  y  -  £■  J]  - 

k  X  (3.11) 

a{x,y)-i{x-/3^,y-sj-  \n{d^  (x,  y )] 


The  Lagrange  constraint  forces  consistency  in  the  original  image  estimation,  z(x,y) , 


iteration  by  iteration,  where  the  Lagrange  can  be  referred  to  as  ; 
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(3.12) 


L{a,  0  =  X  Z  Z  ■  ln[«(x,  y)-i{x-  P^,y-£^)\- 

k  X  y 

a{x,y)-i{x- p,^,y- £^)- XY,'Y4^{x,y) 


X  y 


Now  the  derivative  of  the  log-likelihood  deseribed  in  equation  (3.12)  is  taken  with 
respeet  to  the  gain  nonuniformity  for  one  pixel.  Note  that  the  derivative  of  the  last  term 
in  equation  (3.12)  is  zero  beeause  it  does  not  have  gain: 


d{L{a,i))  _  d 
e(«(Xo,yo))  5(«(Xo,yo)) 

5 


Z  Z  Z  -f )  ■  ln[«(x,  y)-i{x-p„y-£, )] 


\  k  X  y 


a(«(Xo,yo)) 


{a{x,y)-i(x- /3,,y -£,))- 


d 


a(«(Xo,yo)) 


J 

f  \ 

^ZZ'(^’k) 

V  ^  r  J 


(3.13) 


Set  the  result  of  the  derivative  in  equation  (3.13)  equal  to  zero  and  solve  for  the 
gain  at  pixel,  (Xq,  Jo): 


dk(Xo,yo) 


-i{x^-P„y^-£^) 


(3.14) 


Distribute  the  sum  over  k: 


0  = 


Z  ^ki^o’yo) 


(3.15) 


used: 


A  simple  example  with  k  summing  from  1  to  2  demonstrates  the  mathematies 


Q  ^  di(Xo,yo)  ^  d^{x^,y^)  _  -/3^,y^-  s,)-  i{x^  -  l3^,yQ-  s^)  (3.16) 


«(Xo,yo)  «(^o’ko) 

Pull  gain  nonuniformity,  a  ,  out  of  the  sum  over  d  below: 
1 


0  = 


«(^0’ko) 


[d,{x^,y^)  +  d^{x^,y^))-i{x^- p„y^-£^)-i{x^- p^,y^-£^)  (3.17) 
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Move  the  image,  i,  to  the  left  side  of  the  equation: 


'■(^0  -  -  ^i) + '(^0  ^2) 


1 


«(^0>7o) 


(Ji(Xo,Jo)  +  A(^o>7o)) 


(3.18) 


Solve  for  the  gain,  alpha: 


'K-ApJo-^i)  +  'K-A2.Jo-^2) 


(3.19) 


Translate  back  having  sums  over  k  in  the  numerator  and  denominator  to  show  the 
final  result  for  the  gain  nonuniformities: 


«(^o>To)  = 


k _ 


(3.20) 


Now,  substitute  equation  (3.20)  into  equation  (3.12)  to  create  a  log  likelihood 
equation  based  solely  on  the  image  as  shown: 


A0  =  XZZ  d^{x,y)-\n 

k  X  y 


Z  dk^{x,y)-i{x-/3,,y-s,) 

Jh _ 


Z  d,^{x,y)-iix-/3„y-s,) 

"???  “  I  i(^-A,.y-sO 


(3.21) 


Now  differentiate  equation  (3.21)  with  respect  to  the  image  and  set  to  zero: 


d{L{i))  _  dir‘  Term)  ^  8(2’“'  Term)  ^  5(3"''  Term) 


d{i{Y^,ri^))  5(/(yo,?7o)) 


(3.22) 


In  order  to  simplify  the  derivative  operation  the  derivative  of  each  term  is 
demonstrated  individually,  where  the  first  term  using  the  notation  described  above  can  be 
expressed  as: 
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5(r'  Term) 
S{Kro,rio)) 


d 


k  X  y 


S  dk^{x,y)-i{x-/3„y-£,) 
_y _ 

l‘2 


(3.23) 


Where  a  simplifieation  of  the  first  term  results  as  shown  below  using  the 
properties  of  the  natural  log: 


S{i{ro,rio))  k 


dk{x,y)- 


In 


S  dk^{x,y) 


-In 


+  \n[i{x-p„y-£,)] 


(3.24) 


Where  the  derivative  with  respeet  to  the  image  of  the  first  part  inside  the  braekets 
of  equation  (3.24)  above  is  equal  to  zero  beeause  there  is  no  image,  and  simply  a 
derivative  of  a  eonstant: 


diK/o^r/o))  k 


ZZZ  d,(x,y)\n 


=  0 


(3.25) 


The  derivative  of  the  third  part  inside  the  braekets  of  equation  (3.24)  above  is  the 
next  most  ehallenging  and  ean  be  deseribed  as  shown  below  in  equation  (3.26): 

T77-^ — ^ZZZ  d,ix,y)-ln[iix-/3„y-£,)]  = 

S(z(r„.n„  n  V  ,  „ 

(3.26) 


yyy  — ^ 

k  X  ^  S(i(ro,%)) 


<ik{x,y)-\vi[i{x-  P^,y  -  £^)] 


And  ean  be  further  simplified  using  the  ehain  rule: 
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(3.27) 


d,{x,y) 


8 


ZEE  - 

,  ,  j,  i{x-j3„y-s,)d{i{ro,rio)) 


[i{x-/3„y-s,)] 


Where  the  rightmost  statement  in  equation  (3.27)  can  be  calculated  as  follows: 


S(Kro,%)) 


[i{x-p^,y-s^)\  = 


S'EE  ^ir^'n)-S{x-r-(i„y-ri-s^) 


r  1 


5-^(/o.^o) 


(3.28) 


The  derivative  taken  of  the  image  at  any  other  point  other  than  i(7o?7o) 


and  the  single  remaining  point  is: 

50'(ro.7o))  s-i(ro,?io) 


Resulting  in  just  the  shifted  delta  function: 

a 


(3.29) 


Kx-/3k,y-£k)]  =  S{x-ro-/3„y-rio-s,) 


(3.30) 


Making  the  result  of  the  third  part  in  the  brackets  of  equation  (3.24): 
— ^EEE  d,ix,y)-ln[iix-/3„y-s,)]  = 

Siiiro,rio))kx  y 


■y  "y  "y 

k  X  ^  i{x-p^,y-s,) 


(3.31) 


■3{x-y^-p^,y-ri^-s,) 


The  second  part  in  the  brackets  of  equation  (3.24)  can  be  described  as: 


d^{x,y)-\n 


siKn^do))  k  X  >. 

EEE(“i)'^-t(^’k) 


ki 

8 


k  X  y 


-In 


(3.32) 


Where  the  second  part  in  brackets  of  equation  (3.24)  above  is  simplified  to  a  very 
similar  form  found  in  the  third  part  in  brackets  of  equation  (3.24),  having  a  sum  of  the 
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images  over  k  inside  the  natural  log  instead  of  just  one  shifted  image.  Using  this 
similarity  and  the  process  described  above  in  equations  (3.26)-(3.31),  the  final  answer  for 
the  second  part  in  brackets  of  equation  (3.24)  can  be  found  as; 


5(f(ro.^o)) 


*2 


(-Dili 


k  X  y 


ki _ 

ki 


ki 


(3.33) 


Putting  equation  (3.25)  -  (3.33)  together,  the  final  derivative  of  the  first  term  of 
equation  (3.22)  is: 


a(r^  Term) 

S(i(ro’%)) 


zzz 

k  X  y 


dk{x,y) 

Kx-/3,,y-s,) 


■5{x-Y,-p„y-ri^-s,) 


-ZEE 


Z  <^1. 

_ ^ _ 

ki 


•X  d{x-y^-p,^,y-r],-s,^) 

ki 


(3.34) 


The  derivative  of  the  second  term  using  the  notation  described  above  in  equation 


(3.21)  and  equation  (3.22)  is  expressed  as: 


a(2"°'  Term) 


d 

S(i(Yo,%)) 


■(-i)ZZZ 

k  X  y 


X  dk^{x,y)-i{x- p^,y-s^) 

^3 _ 

ki 


(3.35) 


Where  the  shifted  images  can  be  expressed  using  equation  (3.7): 


5(2"'' Term)  _  (-l)-a 

5(f(ro2  7o)) 


ZEE 

k  X  y 


X  dk^{x,y)-Y,Yj  Kr^'n)-S{x-Y-Pk^y-ri-s^) 

Jh _ r  1 _ 

ZEE  i{Y,r/)-S{x-Y-j3,^,y-r/-£,^) 

ki  r  V 


(3.36) 
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The  derivative  with  respect  to  ii/o,riQ)  for  any  image  pixel  outside  is 


zero  making  all  terms  not  at  location  go  to  zero; 


a(2"‘^Term)_  (-1)-^ 

SiKro,rio))  QiKMj) 


III 

k  X  y 


Z  t)  •  '(ro .  ^o)  •  -  7o  -  A .  T  -  '/o  - 

_^3 _ 

h 


(3.37) 


And  can  be  simplified  as: 


djT^  Term) 


(-!)• 


dAx,y)- 


•1 


(3.38) 


Giving  the  simple  result  for  the  second  term  of  zero; 

a(2”^  Term)  _  ^ 
S{i{ro,rio)) 


(3.39) 


The  derivative  of  the  third  term  is  also  rather  simple,  as  there  is  only  a  A  with  the 
sum  of  the  image; 


3(3’'“  Term) 
S(Kro,rio)) 


3 

S(Kn,rio)) 


•(-i)a-ZZ 


(3.40) 


Where  the  derivative  of  sum  over  x  and  y  with  the  image  goes  to  1  and  leaves  the 
constant  A ; 


3(3- Term)  _ ^ 

Giving  the  final  result  for  the  third  term  of  equation  (3.22): 


(3.41) 
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(3.42) 


a(3-Term)^  ^ 

mr,,!,)) 


Putting  the  first,  second  and  third  terms  in  equation  (3.34),  equation  (3.39)  and 
equation  (3.42)  together  to  give  the  final  result; 


d{L{i)) 

SiKn^rio)) 


dkix,y) 

Kx-Pk^y-Sk) 


■s{x-n-pk^y-%-^k) 


S  dkSx,y) 

^2 _ y 

S  Kx-p^,y-sP)  X 

ki 


(3.43) 

5{x-Y^-p„y-ri^-s,)-X 


Using  equation  (3.43)  the  derivative  of  the  log  likelihood  with  respect  to  the 
image  is  taken  and  set  equal  to  zero  maximizing  for  the  image  as  described  above  and 
giving  rise  to  an  iterative  solution; 


k  X  y 


dkix,y) 

Kx-p^,y-sP) 


■dix-Yo-pk^y-rio-^k) 


-ZEE 

k  X  y 


S  dkPx,y) 

^2 _ y 

^2 


5{x-Y^-pk,y-ri^-Sk)-^ 


(3.44) 


Where  the  iterative  solution  for  the  image  using  equation  (3.44)  is; 


hidiMP 


ZEE  -<!,)->■ 

k  X  y  ^\X  Pk^y  ^k) _ 

X  dkSx,y) 

XXX  - o - f'X  d{x-Y^-Pk,y-ri^-Sk) 

k  X  y  2^  i{x-p„y-sP)  k. 


(3.45) 


This  iterative  solution  is  set  up  as  shown  because  the  elements  in  equation  (3.44) 
are  a  positive  gradient,  negative  gradient,  and  the  Lagrange,  and  forming  an  iterative 


-  18- 


solution  that  is  multiplicative  verses  additive  will  ensure  that  the  image  never  goes 
negative,  as  demonstrated  in  Riehardson  [7].  This  is  important  beeause  the  image  is  in 
units  of  photons.  Making  the  multiplier  a  ratio  of  the  positive  gradient  over  the  negative 
gradient  moves  the  old  image  eloser  to  the  ML  estimate  iteration  by  iteration. 

In  order  to  simplify  the  expression  the  equation  ean  be  broken  down  into  what 
will  be  referred  to  in  this  document  as  the  numerator  and  denominator; 


Numerator  =  =  III 

k  X  y 


dk{x,y) 


P^,y-71^-s^)  (3.46) 


Denominator  =  den(yi^,rh^)  = 


III 


I 

ki _ ■y 

S  I3^,y-s^)  X 

h 


5{x-r^-P,,y-ri^-s,) 


(3.47) 


With  these  expressions  the  equation  can  be  simplified  as; 


numiro,ilo)  - 
den(ro,7lo) 


(3.48) 


The  Lagrange  multiplier  needs  to  be  adjusted  for  each  iteration  and  can  be  solved 
for  as  demonstrated  below.  As  shown,  it  is  simplified  in  terms  of  the  numerator  and 
denominator  described  by  equation  (3.46)  and  equation  (3.47); 


LJro,fio)  =  ioidiro,%)- 


numiro,ilo) 

deniro,%) 


ioidiMo) 


A 

deniro,flo) 


(3.49) 


Where  the  resulting  solution  for  the  Lagrange  multiplier  is; 

hidiro,rio)-num{ro,%)  yyy  d,(x,y) 

den(ro,%) _ k.y  ^ 

( ioid(ro,%)  ^ 
yden(ro,%)  y 


(3.50) 
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This  notation  is  especially  helpful  for  coding  practices  as  the  numerator  and  the 
denominator  need  to  be  calculated  for  each  iteration,  then  the  value  for  X  is  found  with 
the  numerator  and  denominator,  and  lastly  the  same  numerator  and  denominator  are  used 
with  the  new/^  in  order  to  find  • 

As  the  iteration  starts,  an  initial  starting  value  for  4z<i(7o’^o)  niust  be  found.  The 
algorithm  will  not  run  to  completion  with  a  continually  decreasing  MSB  if  zeros  or  ones 
are  used  as  starting  place  for  • 

The  first  estimate  for  4z^/(7o’^o)  the  average  result  of  the  frames  with  the  shifts 
included  and  can  be  described  as; 

y)  =  - - - -  (3-51) 


As  a  reminder,  from  equations  (3.11)  through  (3.43)  it  can  be  seen  that  the  log- 
likelihood  with  respect  to  the  gain  nonuniformities  is  found  by  taking  the  derivative  of 
the  log-likelihood  with  respect  to  the  nonuniformities,  setting  it  equal  to  zero,  and  then 
rearranging  the  result  to  find  the  nonuniformities  with  the  result,  as  shown  above; 


«(/o>^o) 


_ k _ 

k 


(3.20) 


This  description  of  the  nonuniformities  can  be  used  to  describe  the  array 
characteristics  as  well  as  used  to  derive  the  stopping  criteria  during  the  iteration  process 
to  find  the  nonuniformities. 
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Requirements  of  Algorithm 


This  algorithm  requires  a  data  eube  with  shifts  that  ean  be  determined  from  the 
scene  data,  or  scene  data  with  the  shifts  in  x  and  j  are  already  provided.  In  the  event 
shifts  are  not  known,  or  cannot  be  found,  this  algorithm  is  not  applicable.  Knowing  this, 
methods  for  shift  estimation  are  used,  however  shift  estimation  is  not  discussed  in  great 
detail  during  the  derivation  of  this  algorithm. 

For  the  purposes  of  testing  this  algorithm  with  measured  data,  a  shift  estimation 
technique  using  two  frames  is  used,  and  can  be  described  generally  by  the  correlation  in 
both  dimensions  between  a  reference  image  and  the  shifted  image,  with  the  resulting 
estimated  shift  being  where  the  correlation  between  frames  was  at  its  peak  in  each 
dimension.  A  reference  for  shift  estimation  can  be  found  in  Robinson  and  Milanfar  for 
additional  information  as  well  as  other  shift  estimation  techniques  [8]. 

In  addition,  the  implementation  requires  that  the  magnitude  of  the  shifts  be  greater 
than  one-half  of  a  pixel  as  the  algorithm  was  derived  solely  with  whole  pixel  shifts  in 
mind,  and  a  pixel  shift  of  magnitudes  one -half  or  less  will  not  come  to  a  desirable 
solution  for  the  original  image.  This  requirement  arises  from  the  need  for  each  pixel  to 
see  a  different  part  of  the  scene  with  the  same  nonuniformity.  An  example  of  the 
dependence  on  accurate  shifts  of  greater  than  one -half  of  a  pixel  in  magnitude  is  shown  in 
the  results  section  and  shows  that  ideally,  whole  shifts  and  accurate  estimation  of  those 
shifts  would  result  in  the  best  performance. 
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It  should  be  noted  that  empty  scenes  (and  the  corresponding  sections  of  the  FPA) 
will  only  contain  bias  data  and  will  not  produce  favorable  results  in  terms  of 
nonuniformity  estimation  and  correction  without  additional  information  from  those  pixels 
in  the  FPA. 

In  terms  of  scene  content,  the  data  cube  to  be  used  should  also  have  image  data  in 
terms  of  received  photons,  and  if  not,  the  data  should  be  converted  to  the  number  of 
received  photons  in  order  to  estimate  nonuniformities,  and  that  the  data  cube  not  contain 
zeros  or  negative  numbers  as  the  estimate  makes  the  assumption  that  there  is  always 
some  bias  for  each  detector. 

It  should  also  be  noted  that  this  algorithm  corrects  for  dead  pixels  (pixels  that 
return  a  near-zero  response  regardless  of  scene  input)  and  naturally  removes  these  dead 
pixels,  assuming  that  the  shifted  data  provides  enough  information  to  correct  for  the  dead 
pixel.  For  example,  if  a  3x3  area  of  pixels  is  dead  and  the  available  image  data  only  has 
shifts  of  up  to  I  pixel  in  any  direction,  the  dead  pixel  area  will  not  be  completely 
corrected  with  this  algorithm. 

Again,  this  algorithm  implementation  is  derived  specifically  for  integer  shifts.  If 
a  data  cube  possesses  subpixel  shifts  the  shift  data  used  in  the  algorithm  should  be 
rounded  to  the  nearest  whole  shift  as  mentioned  above.  This  gives  rise  to  the  requirement 
above  for  the  magnitude  of  shifts  to  be  greater  than  a  half-pixel  in  each  direction  in  order 
to  estimate  the  nonuniformities  in  addition  to  having  the  input  shifts  as  whole  pixels. 
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Derivation  of  Stopping  Criteria 


During  calibration  of  real  LIDAR  data  the  original  seene  is  not  known,  therefore 
there  must  be  some  metrie  for  deciding  when  to  stop  the  iteration  proeess.  This  metrie 
must  use  data  that  is  available  with  a  realistie  seenario.  In  LIDAR,  the  bias  before  signal 
return  is  reeorded  and  the  seene  imagery  itself  with  shifts  must  be  used  to  stop  the 
iteration  to  find  a  satisfaetory  result. 

Assuming  Poisson  noise  in  the  original  scene,  the  mean  of  the  Poisson  noise  is 
equal  to  the  varianee  of  the  same. 

The  average  estimated  noise  varianee  can  be  determined  by; 

,  (x,  y )  -  (x,  y )  •  (x,  y ) )" 

Noise  Variance  =  — — ^ — -  (3.52) 

K-X-Y 

And  the  mean  ean  be  deseribed  by  the  following; 

Mean=  ^  ^  ^ -  (3.53) 

K-X-Y 

If  the  average  variance  is  less  than  the  average  mean  the  iteration  is  terminated 
beeause  the  average  varianee  shows  that  an  improvement  have  been  made  to  the  image, 
giving  rise  to  the  stopping  eriteria; 

k  X  y  k  X  y 
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Where  (x„^„{x,y)  are  the  estimated  nonuniformities,  and  K, 

number  of  frames,  the  length  of  x,  and  the  length  of  y,  respeetively. 
equation  beeomes  true  the  algorithm  terminated. 


X,  and  Y  are  the 
When  the  above 
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IV.  Results  of  Algorithm 


In  order  to  test  the  new  algorithm  a  simple  noiseless  case  was  created  with  a 
simple  ramp  scene  in  one  dimension.  After  the  algorithm  worked  with  the  noiseless  case, 
Poisson  noise  was  added  to  the  one  dimensional  ramp  and  the  algorithm  was  tested  again. 
Afterward,  a  two  dimensional  ramp  case  was  tested  without  Poisson  noise,  and  then  with. 
For  completeness,  another  two  dimensional  test  scene  was  created  with  sinusoids  and  is 
shown  below.  These  cases  were  necessary  in  order  to  test  and  validate  the  model  and  test 
derived  stopping  methods.  In  all  cases  the  gain  nonuniformities  varied  uniformly  from 
0.9  to  1.1.  A  brief  study  of  the  algorithm’s  dependence  on  accurate  shift  estimates  is  also 
included  at  the  end  of  this  chapter. 


ID  No  Noise  Performance 

A  one  dimensional  noiseless  case  was  the  first  test  applied  to  the  derived 
algorithm.  This  case  was  used  as  a  test  method  and  has  aided  significantly  in  algorithm 
development  and  testing.  This  one  dimensional  case  contained  a  scene  with  a  simple 
intensity  ramp  and  low  values  on  the  outer  edges.  The  maximum  value  for  the  ramp  was 
300,000  units,  based  on  a  value  that  would  approximate  a  maximum  well  depth  (Well 
depth  being  the  maximum  number  of  photons  before  a  readout  is  necessary.)  for  photons 
in  a  sensor.  The  final  image  fed  to  the  algorithm  consisted  of  the  original  scene  data 
multiplied  by  the  nonuniformities  as  shown  below; 

dk  {x,  y)  =  i{x-  ijvix,  y)  +  n{x,  y)  (4.1) 
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Where  d,^{x)  represents  the  image  with  gain  nonuniformities, 
represents  the  original  scene  shifted  by  the  A:*  shift  in  x,  and/7L'(x)  represents  the  gain 
nonuniformities  in  the  focal  plane  array. 

The  performance  of  the  algorithm  under  the  one  dimensional  noiseless  case  is 
extremely  satisfactory.  With  only  two  image  frames,  and  a  one  pixel  shift,  a  very  close 
approximation  to  the  original  scene  with  nonuniformities  was  found,  showing  the  noisy 
data,  or  data  with  only  nonuniformities  present.  The  results  can  be  seen  in  figure  1 
below; 


Figure  2  -  Results  of  corrected  image  compared  to  original  with  nonuniformities. 
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It  can  be  seen  in  figure  2  above  that  the  corrected  one  dimensional  image  is  a 
much  closer  match  to  the  true  image  than  the  noisy  data  (data  with  nonuniformities,  no 
Poisson  noise).  The  shift  between  the  two  frames  used  in  this  example  is  simply  1  and 
the  result  in  figure  2  and  3  resulted  after  50  iterations. 

As  the  algorithm  is  iterated  the  MSB  (Mean  Squared  Error)  is  reduced  iteration  by 
iteration  and  eventually  bottoms  out,  as  shown  in  figures  3  and  5; 


Figure  3  -  Running  MSB  (Mean  Square  Error)  iteration  by  iteration. 

Figure  3  above  is  a  Bog-Bog  figure  corresponding  to  the  MSB  for  figure  2,  it 
appears  from  a  limited  number  of  iterations  the  MSB  appears  to  go  down  as  the  iteration 
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number  increases,  but  after  many  iterations  the  MSB  does  reach  a  minima,  as  shown  in 
figure  5. 

With  increased  iterations  the  performance  of  the  algorithm  continues  to  improve 
the  overall  image  and  the  MSB  continues  to  be  reduced,  as  shown  in  figure  4; 


figure  4  -  One  dimensional  test  case  with  many  iterations  and  no  Poisson  noise. 


It  can  be  seen  in  figure  4  that  with  even  more  iterations  the  resulting  image 
quality  is  improved  even  more.  This  is  especially  evident  when  figure  4  is  compared  to 
the  error  seen  in  figure  2  in  the  upper  right  corner. 
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As  the  algorithm  is  iterated  to  an  even  higher  number  the  MSB  is  reduced  even 
further  and  eventually  bottoms  out.  This  can  be  seen  in  the  following  figure; 


fl  1  5  3  4 
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Iteration  Number 

Figure  5  -  Resulting  MSB  from  no  Poisson  noise  case  with  many  iterations. 

It  can  be  seen  by  the  log-log  plot  in  figure  5  that  the  corresponding  to  the  MSB 
that  results  after  producing  figure  4  that  the  MSB  is  reaches  a  minima  after  roughly  5000 
iterations  and  runs  to  roughly  10,000  iterations  as  shown  in  figure  4. 

While  it  may  not  be  evident  in  figure  3  that  the  MSB  contains  a  minima,  with  only 
quantization  noise  and  nonuniformities  the  MSB  does  reach  a  minima  after  roughly  5000 
iterations  as  shown  in  figure  5.  Note  that  this  would  be  the  optimum  stopping  point  for 
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the  algorithm.  Figures  with  Poisson  noise  hereafter  will  demonstrate  the  stopping 
criteria. 


ID  With  Poisson  Noise  Performance 

After  the  testing  without  noise  was  completed  the  algorithm  was  tested  with 
Poisson  noise.  The  image  with  Poisson  noise  was  then  processed  by  the  algorithm  and 
the  results  can  be  seen  below  in  figure  6; 


Figure  6  -  Resulting  image  with  Poisson  noise  after  5000  iterations. 
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It  can  be  seen  in  figure  6  that  with  5000  iterations  the  resulting  image  quality 
improved,  even  with  the  presence  of  Poisson  noise,  though  not  nearly  as  well  as  in  figure 
4. 

The  resulting  MSB  is  also  decreased  with  the  presence  of  Poisson  noise.  It  should 
be  noted  the  presence  of  Poisson  noise  prevents  the  same  algorithm  from  reducing  the 
MSB  in  the  same  number  of  iterations,  where  the  Poisson  noise  causes  the  algorithm  to 
end  earlier  and  results  in  less  correction  than  the  case  without  Poisson  noise  present; 
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Bigure  7  -  The  MSB  resulting  from  a  test  case  with  Poisson  noise. 
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In  the  log-log  plot  on  figure  7  it  can  be  seen  with  that  the  algorithm  does  not 
perform  as  well  with  Poisson  noise  than  without  it.  It  can  also  be  seen  that  the  MSB 
bottoms  out  after  many  iterations,  with  the  stopping  criteria  ending  at  4000  iterations 
while  the  optimum  stopping  point  would  be  at  3000  iterations. 

Figure  7  reveals  characteristics  of  the  MSB  verses  iteration  number,  showing  that 
the  image  reaches  a  best  case  with  noise  and  then  if  the  iteration  continues  the  resulting 
image  quality  is  reduced.  If  allowed  to  continue  the  MSB  will  continue  to  rise. 


2D  No  Noise  Performance 

For  the  two  dimensional  case  without  Poisson  noise  a  two  dimensional  ramp 
similar  to  the  one  dimensional  case  with  a  maximum  value  of  300,000  units  at  the 
maximum  corner  and  near-zero  at  the  minimum  corner.  Just  as  with  the  no  Poisson  noise 
one  dimensional  case,  the  two  dimensional  case  produces  very  favorable  results  as  seen 
by  a  slice  of  the  image  through  the  50**'  row: 
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Pixel  Return  (Relative) 


Figure  8  -  2D  ramp  without  Poisson  noise. 

Figure  8  shows  a  slice  through  the  two  dimensional  ramp  it  can  be  seen  in  the 
corrected  one  dimensional  image  that  the  dots  representing  the  estimated  image  are  a 
much  closer  matching  to  the  true  image  than  the  noisy  data  with  nonuniformities. 

The  resulting  MSB  is  also  very  favorable  as  evidenced  by  the  following  plot; 
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Figure  9  -  The  resulting  MSB  for  the  2  dimensional  ramp  without  Poisson  noise. 

In  figure  9  the  MSB  appears  to  be  decreasing  iteration  by  iteration.  As  with  the 
one  dimensional  case  the  MSB  does  bottom  out  after  many  thousands  of  iterations, 
although  not  obvious  from  this  figure. 

The  nonuniformities  in  the  images  themselves  and  the  corrected  image  using  the 
algorithm  can  be  seen  by  the  comparative  illustration  below; 
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With  NUs  NUs  Corrected 


Figure  10  -  An  illustration  with  the  2D  test  ramp  with  and  without  nonuniformities. 

Figure  10  shows  a  eomparative  illustration  of  the  test  ramp  with  nonuniformities 
showing  the  image  with  nonuniformities  on  the  left  and  the  eorreeted  image  on  the  right. 
Note  that  there  is  no  Poisson  noise  in  this  case,  giving  a  very  favorable  result. 


2D  With  Noise  Performance 

For  the  two  dimensional  noise  case  a  two  dimensional  ramp  similar  to  the  one 
dimensional  case  with  a  maximum  value  of  300,000  units  at  the  maximum  corner  and 
near-zero  at  the  minimum  corner  as  mentioned  above.  Using  this  same  image,  Poisson 
noise  is  included  in  each  frame  of  the  image  series  and  the  resulting  frames  are  processed 
by  the  algorithm.  Just  as  with  the  one  dimensional  case  with  Poisson  noise,  the  two 
dimensional  case  produces  favorable  results  as  seen  by  a  slice  of  the  image  through  the 
50*  row  of  the  image; 
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Figure  1 1  -  2D  ramp  with  Poisson  noise. 


The  image  slice  in  figure  1 1  above  demonstrates  a  portion  of  the  corrected  two 
dimensional  image  with  Poisson  noise  present.  Note  the  lower  performance  compared  to 
the  case  without  Poisson  noise. 

The  resulting  MSB  is  also  very  favorable  as  evidenced  by  the  following  plot  with 
the  stopping  criteria  ending  slightly  prematurely  but  very  close  to  the  minimum  MSB; 
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Figure  12  -  The  MSB  resulting  from  the  2D  test  ramp  with  Poisson  noise. 

The  resulting  log-log  MSB  in  figure  12  above  shows  excellent  performance 
iteration  by  iteration  until  it  reaches  a  minima  after  roughly  1000  iterations  while  the 
stopping  criteria  stops  the  iteration  after  roughly  200  iterations  with  this  relatively  low 
level  of  Poisson  noise. 

The  nonuniformities  in  the  images  themselves  and  the  corrected  image  using  the 
algorithm  can  be  seen  by  the  comparative  illustration  below; 
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With  NUs  Only  With  NUs  and  Noise  NUs  Corrected  with  Noise 

Figure  13  -  An  illustration  showing  images  with  NUs,  NUs  and  Poisson  noise,  and  the 

final  eorreeted  image. 

The  eomparative  illustration  in  figure  13  above  shows  the  nonuniformities 
ineluding  two  dead  pixels,  the  nonuniformities  with  Poisson  noise,  and  the  resulting 
eorreeted  image  showing  the  dead  pixels  fdtered  out.  Both  dead  pixels  are  centered 
horizontally,  one  near  the  top  and  the  second  near  the  bottom. 

Additional  2D  Case  With  Noise 

In  order  to  demonstrate  that  the  algorithm  works  on  more  than  a  simple  ramp, 
another  case  was  prepared  to  demonstrate  the  algorithm’s  effectiveness.  The  figure 
below  shows  the  performance  of  this  test  image  with  a  two  dimensional  sinusoid; 
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Figure  14  -  2D  sinusoidal  test  case  with  Poisson  noise. 

It  can  be  seen  in  figure  14  above  that  the  image  slice  of  a  two  dimensional 
sinusoidal  pattern  that  the  resulting  image  is  corrected. 

It  should  also  be  noted  that  the  resulting  MSB  shows  that  the  algorithm  bottoms 
out  quickly,  and  slowly  increases  iteration  by  iteration  after  the  minima  is  reached; 
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Figure  15  -  MSB  that  results  from  2D  sinusoidal  test  case  with  Poisson  noise. 

The  resulting  MSB  in  figure  15  bottoms  out  quickly,  and  slowly  rises  over  many 
iterations  where  the  estimated  stopping  point  stops  at  roughly  1950  iterations  with  only  a 
minimal  increase  in  the  MSB  still  providing  an  excellent  correction. 

The  resulting  images  from  this  two  dimensional  sinusoidal  test  case  are  shown 
below  where  the  location  of  two  dead  pixels  are  easily  seen  in  Bigure  16; 
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Figure  16  -  An  illustration  showing  the  2D  sinusoidal  test  ease  with  NUs,  NUs  and 
Poisson  noise,  and  the  corrected  image. 


It  can  be  seen  in  figure  16  above  that  the  two  dimensional  sinusoidal  image,  even 
with  two  dead  pixels,  is  corrected  and  the  dead  pixels  are  removed. 

Overall,  the  performance  of  the  algorithm  on  simulated  images  was  excellent, 
with  the  MSB  reaching  its  low  point  very  quickly,  and  the  dead  pixels  being  removed 
from  the  original  image. 


Algorithm  Dependence  on  Accurate  Shift  Estimates 

The  algorithm  was  derived  using  whole  shifts,  and  as  a  result  is  dependant  on 
accurate  shift  inputs  and  their  related  input  scenes.  In  order  to  test  the  dependence  on 
accurate  shifts,  a  set  of  2D  test  image  pairs  with  sub-pixel  shifts  between  0.6  and  1.4 
pixels  in  x  andy  were  generated.  Each  pair  of  images  was  used  with  an  input  of  a  shift  in 
X  and  y  and  the  MSE  results  were  compared.  The  results  are  shown  in  the  figures 
following; 
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Comparison  of  subpixel  shifts  as  an  input  to  algorithm 
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Figure  17  -  A  demonstration  of  algorithm  sensitivity  to  accurate  shift  inputs. 

Figure  17  above  demonstrates  sets  of  image  pairs  with  real  sub-pixel  shifts  in  x 
andy  ranging  from  0.6  to  1.4  pixels.  The  minimum  MSB  achieved  in  each  scenario  is 
plotted  showing  that  the  minimum  MSB  is  reached  when  the  shift  estimate  fed  to  the 
algorithm  is  closest  to  the  accurate  value  of  the  shift.  Notice  the  higher  error  on  shifts  of 
1.1  to  1.4  because  of  the  nature  of  the  ramp  image  used  to  test  the  scenario. 

Notice  that  the  minimum  MSB  is  reached  when  the  shift  data  fed  to  the  algorithm 
is  most  accurate.  This  is  expected  as  the  algorithm  was  derived  to  use  whole  shifts  in 
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order  to  optimize  for  speed.  The  error  introduced  by  uncertainty  in  the  shifts  is  highly 
dependant  on  the  original  image  data.  This  can  be  seen  by  the  higher  error  on  the  right 
side  of  the  graph  on  the  figure  above  due  to  the  original  image  being  a  ramp  that 
increases  to  the  right  and  down,  and  a  resulting  error  in  the  shift  down  and  to  the  right 
increases  the  relative  error  between  the  scenarios. 

Additionally,  the  running  MSB  was  compared  in  each  scenario  in  a  log-log  plot, 
showing  that  the  MSB  performs  well  iteration  by  iteration  when  the  shift  estimate  fed  to 
the  algorithm  is  most  accurate; 


Bigure  18  -  Running  MSB  by  iteration  of  inaccurate  shift  inputs. 


-43  - 


Figure  18  demonstrates  the  running  MSB  iteration  by  iteration  with  sets  of  image 
pairs  with  sub-pixel  shifts  in  x  and  y  ranging  from  0.6  to  1.4  pixels.  The  minimum 
running  MSB  is  achieved  when  the  shift  estimate  fed  to  the  algorithm  is  closest  to  the 
accurate  value  of  the  shift.  Notice  the  higher  error  on  shifts  of  1 . 1  to  1.4  shown  in 
magenta  because  of  the  nature  of  the  ramp  image  used  to  test  the  scenario. 

The  log-log  plot  in  figure  1 8  also  demonstrates  the  eventual  bottoming  out  of  the 
MSB,  especially  when  incorrect  shifts  are  inputted  into  the  algorithm.  Also  note  that 
while  it  is  not  especially  visible  in  figure  18,  the  1.0  pixel  shift  case  does  bottom  out  after 
many  thousands  of  iterations. 
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V.  Conclusion 


Random  gain  nonuniformities  in  LIDAR  arrays  are  a  problem  because  of  the 
nature  of  these  arrays.  Because  of  the  nature  of  LIDAR  arrays,  having  bias  data  for  each 
collection,  the  gain  can  be  corrected  with  the  iterative  solution  described  in  this  thesis.  In 
order  for  this  algorithm  to  operate  as  implemented,  there  must  be  at  least  a  full  pixel  of 
shifts  in  the  horizontal  and/or  vertical  directions  between  two  image  frames  used  as  input, 
and  there  must  be  an  accurate  estimate  of  shifts,  preferably  as  close  to  a  whole  pixel  as 
possible. 

Having  the  ability  to  correct  these  nonuniformities  with  a  low  number  of  frames 
makes  the  algorithm  relatively  fast  in  terms  of  computation  time.  This  also  makes  this 
algorithm  useful  to  the  LIDAR  community,  improving  resulting  images. 

In  test  scenarios,  the  algorithm  performs  quite  well,  especially  in  low-noise  cases, 
and  iterates  to  a  solution  that  is  more  favorable  than  the  original  image,  especially  when 
the  algorithm  is  stopped  at  the  appropriate  iteration.  In  the  event  that  the  algorithm  is  not 
stopped  at  the  appropriate  iteration  the  resulting  image  could  be  less  favorable  than 
desired  showing  little  improvement. 

It  is  also  important  to  note  that  all  the  test  cases  were  performed  knowing  the 
exact  shift  between  image  pairs  tested,  and  those  shifts  inx  andy  were  whole  pixel  shifts, 
with  the  exception  of  the  comparison  of  sub-pixel  shift  performance  where  the  algorithm 
dependence  on  accurate  shift  estimation  is  demonstrated. 

It  is  recommended  that  for  future  work  that  the  magnitude  of  error  resulting  from 
sub  pixel  shifts  is  found  for  real  data  and  the  relationship  between  sub  pixel  shifts  and 
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algorithm  performance  is  discovered  for  more  than  just  simulated  scenarios.  It  is  also 
recommended  a  sub-pixel  version  of  the  algorithm  described  in  this  document  be  derived 
and  tested  in  order  to  compare  performance  and  computational  time  difference.  Work 
demonstrating  that  the  algorithm  described  here  could  operate  in  real  time  along  with  an 
operational  LIDAR  system  could  show  the  additional  usefulness  of  this  algorithm. 
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